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ABSTRACT 

In  this  paper,  we  introduce  a  new  millimeter  wave  imaging  modal¬ 
ity  with  extended  depth-of-field  that  provides  diffraction  limited  im¬ 
ages  based  on  a  significant  reduction  in  scan-time.  The  technique 
uses  a  cubic  phase  element  in  the  pupil  of  the  system  and  a  non¬ 
linear  recovery  algorithm  to  produce  images  that  are  insensitive  to 
object  distance.  We  present  experimental  results  that  validate  sys¬ 
tem  performance  and  demonstrate  a  greater  than  four-fold  increase 
in  depth-of-field  with  a  reduction  in  scan-time  by  a  factor  of  at  least 
two. 

Index  Terms —  Computational  imaging,  millimeter  wave  imag¬ 
ing,  extended  depth-of-field,  image  reconstruction,  sparsity. 

1.  INTRODUCTION 

Over  the  past  several  years,  imaging  using  millimeter  wave  (mmW) 
and  terahertz  technology  has  gained  a  lot  of  interest  [1],  [2],  [3].  This 
interest  is,  in  part,  driven  by  the  ability  to  penetrate  poor  weather  and 
other  obscurants  such  as  clothes  and  polymers.  Millimeter  waves  are 
high-frequency  electromagnetic  waves  usually  defined  to  be  in  the 
30  to  300  GHz  range  with  corresponding  wavelengths  between  10 
to  1mm.  Radiation  at  these  these  frequencies  is  non-ionizing  and  is 
safe  to  use  on  people.  Applications  of  this  technology  include  the 
detection  of  concealed  weapons,  explosives  and  contraband.  Fig.  1 
compares  a  visible  image  and  corresponding  94-GHz  image  of  two 
people  with  various  weapons  concealed  under  clothing.  Note  that 
concealed  weapons  are  clearly  detected  in  the  mmW  image. 

Recently,  in  [3],  Mait  et  al.  presented  a  computational  imaging 
method  to  extend  the  depth-of-field  of  a  passive  mmW  imaging  sys¬ 
tem.  The  method  uses  a  cubic  phase  element  in  the  pupil  plane  of 
the  system  to  render  system  operation  relatively  insensitive  to  object 
distance.  The  aberrations  introduced  by  the  cubic  phase  elements  are 
then  removed  by  post-detection  signal  processing.  It  was  shown  that, 
one  can  increase  the  depth-of-field  of  a  94GHz  imager  four  times  its 
normal  depth-of-field  [3]. 

Several  other  systems  have  also  been  developed  and  discussed 
in  the  literature  [1],  [2],  Some  of  them  use  a  single-beam  system 
that  forms  an  image  by  scanning  in  azimuth  and  elevation.  One  of 
the  main  drawbacks  of  mechanical  scanning  is  that  it  significantly 
limits  the  acquisition  speed.  For  instance,  it  takes  about  15  seconds 
for  a  94-GHz  imager  to  scan  a  30°-by-30°  angular  sector.  Real¬ 
time  mmW  imaging  has  also  been  demonstrated  using  an  array  of 
sensors.  However,  these  systems  introduce  higher  complexity  and 
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are  costly.  To  deal  with  this,  compressive  sampling  methods  [4],  [5] 
have  been  applied  to  mmW  imaging  which  reduces  the  number  of 
samples  required  to  form  an  image  [6],  [7],  [8],  [9],  [10],  [11]. 


(a)  (b) 

Fig.  1.  Millimeter  wave  imaging  through  clothing,  (a)  Visible  image 
of  the  scene,  (b)  Image  produced  using  a  94-GHz  imaging  system. 


In  this  paper,  we  describe  a  passive  mmW  imaging  system  with 
extended  depth-of-field  that  can  produce  images  with  reduced  num¬ 
ber  of  samples.  We  show  that  if  the  mmW  image  is  assumed  to  be 
sparse  in  some  transform  domain,  then  one  can  reconstruct  a  good 
estimate  of  the  image  using  this  new  image  formation  algorithm.  Our 
method  relies  on  using  a  far  fewer  number  of  measurements  than  the 
conventional  systems  do  and  can  reduce  the  scan-time  significantly. 

The  organization  of  the  paper  is  as  follows.  Section  2  provides 
background  information  on  passive  mmW  imaging  using  a  94GHz 
system.  The  proposed  undersampling  scheme  is  described  in  Sec¬ 
tion  3.  We  demonstrate  experimental  results  in  Section  4  and  Sec¬ 
tion  5  concludes  the  paper  with  a  brief  summary  and  discussion. 

2.  BACKGROUND 

A  schematic  diagram  of  our  mmW  imaging  system  is  shown  in 
Fig.  2.  It  is  a  94-GHz  Stokes-vector  radiometer  used  for  mmW 
measurements.  It  is  a  single-beam  system  that  produces  images  by 
scanning  in  horizontal  and  vertical  axis.  The  radiometer  has  a  ther¬ 
mal  sensitivity  of  0.3  K  with  a  30-ms  integration  time  and  1-GHz 
bandwidth  per  pixel.  A  Cassegrain  antenna  is  mounted  to  the  front 
of  the  radiometer  receiver  that  has  24”  diameter  primary  parabolic 
reflector  and  a  1.75”  diameter  secondary  hyperbolic  reflector.  The 
position  of  the  hyperbolic  secondary  is  variable. 

One  can  model  the  94-GHz  imaging  system  as  a  linear,  spatially 
incoherent,  quasi-monochromatic  system  [3].  The  intensity  of  the 
detected  image  can  be  represented  as  a  convolution  between  the  in¬ 
tensity  of  the  image  predicted  by  the  geometrical  optics  with  the 


parabolic  primary 


(a)  (b) 


Fig.  2.  94-GHz  imaging  system,  (a)  Image  of  the  scanning  system, 
(b)  Schematic  diagram  with  measured  dimensions. 

system  point  spread  function  [12] 

ii(x,  y)  =  | i(x,  y)\2  =  og(x,  y)  *  *h(x,  y),  (1) 

where  **  represents  the  two-dimensional  convolution  operation. 
The  function  og(x,y)  represents  the  inverted,  magnified  image  of 
the  object  that  a  ray-optics  analysis  of  the  system  predicts. 

The  second  term  in  (1),  h(x,y),  is  the  incoherent  point  spread 
function  (PSF)  that  accounts  for  wave  propagation  through  the  aper¬ 
ture 

h(x-v)-Wr\p{jT-Tf)\  •  (2> 

where  p(x / \f,y / \f)  is  the  coherent  point  spread  function.  The 
function  p(x ,  y)  is  the  inverse  Fourier  transform  of  the  system  pupil 
function  P(u,  v) 

p{x,y)  =  FT_1[P(m,u)]. 


Wv  along  the  respective  axis.  The  constant  7  represents  the  strength 
of  the  cubic  phase.  Fig.  3  shows  the  measured  PSFs  for  conventional 
imaging  and  imaging  with  a  cubic  phase.  Note  that  the  response  of 
the  cubic  phase  system  is  relatively  unchanged,  whereas  the  response 
of  the  conventional  system  changes  considerably. 
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Fig.  3.  Measured  point  spread  functions  for  conventional  imaging 
and  imaging  with  a  cubic  phase.  PSFs  for  conventional  system  at  (a) 
113”,  (b)  146.5”,  and  (c)  180”.  (d)-(f)  PSFs  for  a  system  with  cubic 
phase  at  the  same  distances  for  (a)-(c). 

A  post-detection  signal  processing  step  is  usually  employed  to 
remove  the  artifacts  of  the  aberrations  introduced  by  the  cubic  phase 
element  to  produce  a  well-defined  sharp  response  [13],  [14],  [15]. 
Assuming  a  linear  process 

ip(x,  V )  =  V )  *  *w(x,  y)  (5) 


Without  loss  of  generality,  assume  the  recorded  arrays  are  of  size 
N  x  N.  Then,  eq.  (1)  can  be  described  as 

i  =  Ho9,  (3) 

where  i  and  og  are  N2  x  1  lexicographically  ordered  column  vectors 
representing  the  N  x  N  arrays  ii(x,y)  and  og(x,y),  respectively, 
and  FI  is  the  N2  x  N 2  matrix  that  models  the  incoherent  point  spread 
function  h(x,y).  Displacement,  de,  of  an  object  from  the  nominal 
object  plane  introduces  a  phase  error,  9e  (u,  v),  in  the  pupil  function. 
The  phase  error  increases  the  width  of  a  point  response. 

The  system’s  depth-of-field  {DoF)  is  defined  as  the  distance  in 
object  space  over  which  an  object  can  be  placed  and  still  produce 
an  in-focus  image.  For  a  94  GHz  imager  with  an  aperture  diameter 
D  =  24”  and  object  distance  d0  =  180”,  DoF  «  17.4”  which 
ranges  from  175.2”  to  192.6”  [3], 

The  cubic  phase  element  Pc(u,  v)  is 


PC(U,V) 


exp (jOc{u,  v))rect 


(  u  v  \ 

VWV  Wv)  ' 


(4) 


where 


9c(u,v)  =  (777) 


and  rect  is  the  rectangular  function.  The  phase  function  is  separa¬ 
ble  in  the  u  and  v  spatial  frequencies  and  has  spatial  extent  Wu  and 


one  can  implement  w(x,  y)  as  a  Wiener  filter  in  the  Fourier  space 


W{u,v) 


\Hc(u,v)\2  + 


K—2&n(u,v)  5 


where  Hc(u,v)  is  the  optical  transfer  function  associated  with  the 
cubic  phase  element,  the  parameter  K  is  a  measure  of  the  signal-to- 
noise  ratio,  the  functions  T>£  and  $  at  are  the  expected  power  spec¬ 
tra  of  the  object  and  noise,  respectively.  The  optical  transfer  func¬ 
tion  is  usually  estimated  from  the  experimentally  measured  point 
responses.  One  can  view  the  estimated  ip  as  a  diffraction  limited 
response.  In  the  matrix  notation,  one  can  rewrite  eq.(5)  as 


ip  =  Wi  =  WHo3,  (7) 

where  ip  is  the  N 2  x  1  column  vector  corresponding  to  array  ip(x,  y) 
and  W  is  the  N2  x  N2  convolution  matrix  corresponding  to  the 
Wiener  filter  w(x,y). 


3.  ACCELERATED  IMAGING  WITH  EXTENDED 
DEPTH-OF-FIELD 

Since  our  objective  is  to  form  mmW  images  with  reduced  number  of 
samples,  we  propose  the  following  sampling  strategy.  Our  sensor  is 
a  single-beam  system  that  produces  images  by  scanning  in  azimuth 
and  elevation.  One  can  reduce  the  number  of  samples  by  randomly 


undersampling  in  both  azimuth  and  elevation.  Mathematically,  this 
amounts  to  introducing  a  mask  in  eq.  (1)  as  follows 

iM  =  Mi  =  MHo9,  (8) 

where  i m  is  an  N2  x  1  lexicographically  ordered  column  vector  of 
observations  with  missing  information.  Here,  M  is  the  degradation 
operator  that  removes  the  p  samples  from  the  signal.  This  matrix  is 
of  size  N2  x  N2,  built  by  taking  the  N2  x  N2  identity  matrix  with 
p  zeros  in  the  diagonal  correspond  to  the  discarded  samples. 

To  remove  the  artifacts  of  the  aberrations  introduced  by  the  cu¬ 
bic  phase  element  a  Wiener  filter  can  be  implemented  as  shown  in 
eq.  (5).  Then,  the  observation  model  can  be  written  as 

i0  =  WiM  =  WMi  =  WMHoj.  (9) 

Using  the  relation  in  eq.  (7),  one  can  write  Hos  in  terms  of  the 
diffraction  limited  response,  ip,  as 

Ho9  =  Gip,  (10) 

where,  G  is  the  regularized  inverse  filter  corresponding  to  W.  With 
this,  and  assuming  the  presence  of  additive  noise  77,  one  can  rewrite 
the  observation  model  (9)  as 

io  =  WMGip  +  77,  (11) 

where  r)  denotes  the  N2  x  1  column  vector  corresponding  to  noise, 
77.  We  assume  that  ||t7||2  =  e2. 

Having  observed  i0  and  knowing  the  matrices  W,  M  and  G,  the 
general  problem  is  to  estimate  the  diffraction  limited  response,  ip. 
Assume  that  \v  is  sparse  or  compressible  in  a  basis  or  frame  \P  so  that 
ip  =  vDa  with  ||a||o  =  K  <C  N2,  where  the  do  sparsity  measure 
|.||o  counts  the  number  of  nonzero  elements  in  the  representation. 
The  observation  model  (11)  can  now  be  rewritten  as 

i0  =  WMG$q  +  tj.  (12) 

This  is  a  classic  inverse  problem  whose  solution  can  be  obtained  by 
solving  the  following  optimization  problem 

a  =  argmin  ||  a  ||i  subject  to  ||i0  —  WMGfPa:  II2  <  e.  (13) 

Ot 

Once  the  representation  vector  a  is  estimated,  we  obtain  the  final 
estimate  of  ip  as  ip  =  (Da.  Note  that  the  recovery  of  a  from  eq.  (12) 
will  depend  on  certain  conditions  on  the  sensing  matrix  WMG$ 
and  the  sparsity  of  a  [16]. 

4.  EXPERIMENTAL  RESULTS 

In  this  section,  we  demonstrate  the  performance  and  applicability 
of  our  method  on  simulated  data  with  the  measured  PSFs.  In  these 
experiments,  we  use  an  orthogonal  wavelet  transform  (Daubechies 
4  wavelet)  as  a  sparsifying  transform.  There  has  been  a  number  of 
approaches  suggested  for  solving  optimization  problems  such  as  eq. 
(13).  In  our  approach,  we  employ  a  highly  efficient  algorithm  that  is 
suitable  for  large  scale  applications  known  as  the  Gradient  Projection 
for  Sparse  Reconstruction  (GPSR)  algorithm  [17], 

The  extended  object  used  in  our  experiments  is  represented  in 
Fig.  4(a).  Images  of  an  extended  object  for  conventional  imaging 
system  at  113”,  146”  and  180”  are  shown  in  Fig.  5(a)-(c),  respec¬ 
tively.  Note  that  the  conventional  imaging  system  produces  images 
with  significant  blurring.  In  contrast,  even  without  signal  process¬ 
ing,  the  images  produced  with  cubic  phase  element  retain  more  dis- 
cernable  characteristics  of  the  object  than  the  images  from  the  con¬ 
ventional  system,  as  shown  in  Fig.  5(d)-(f).  It  can  be  seen  from 


Fig.  5(g)-(i)  that  post  processing  removes  the  artifacts  of  aberrations 
introduced  by  the  cubic  phase  element  and  produces  sharp  images. 
Hence,  one  can  extend  the  region  over  which  the  system  generates 
diffraction  limited  images.  In  fact,  in  [3],  it  was  shown  that  the 
DoF  of  a  conventional  94GHz  imaging  system  can  be  extended 
from  17.4”  to  more  than  68”. 


12" 


(a)  (b) 


Fig.  4.  (a):  Representation  of  the  extended  object  used  to  compare 
conventional  and  cubic-phase  imaging,  (b):  Schematic  of  object  il¬ 
lumination. 


Fig.  5.  Images  from  a  conventional  imaging  system  at  (a)  113”,  (b) 
146”  and  (c)  180”.  (d)-(f)  Images  from  a  system  with  cubic  phase 
at  the  same  object  distances  as  for  (a)-(c).  (g)-(f)  Processed  images 
from  a  system  with  cubic  phase  at  the  same  object  distances  as  for 
(a)-(c). 


Fig.  7(a)-(c)  show  the  sparsely  sampled  cubic  phase  data.  Only 
50%  of  the  data  is  sensed.  The  samples  were  discarded  according 
to  a  random  undersampling  pattern  shown  in  Fig.  6(a).  This  corre¬ 
sponds  to  a  reduction  in  scan-time  by  a  factor  of  2.  The  reconstructed 
images  obtained  by  solving  problem  (13)  are  shown  in  Fig.  7(d)-(f). 


The  reconstructions  of  the  extended  object  are  comparable  to  the 
processed  images  from  a  system  with  cubic  phase.  This  can  be  seen 
by  comparing  Fig.  5(g)-(i)  with  Fig.  7(d)-(f). 

In  Fig.  6(b),  shows  the  peak-signal-to-noise-ratio  (PSNR)  curves 
as  we  vary  the  number  of  measurements.  We  see  that  with  increase 
in  number  of  measurements  the  reconstruction  quality  generally  im¬ 
proves.  It  is  interesting  to  note  that  the  PSNR  values  stay  about  the 
same  after  50%  of  the  measurements  are  taken.  Furthermore,  recon¬ 
struction  curves  corresponding  to  all  three  distances  113”,  146”  and 
180”  follow  the  similar  curves,  showing  the  depth  invariance  of  the 
system.  These  figure  shows  that  it  is  indeed  possible  to  extended 
depth-of-field  even  when  cubic  phase  data  is  sparsely  sampled. 


Fig.  6.  (a)  A  random  undersampling  pattern,  (b)  Relative  error  vs. 
number  of  measurement  curves. 


Fig.  7.  Sparsely  sampled  data  from  a  modified  imaging  system  at 
(a)  113”,  (b)  146”,  and  (c)  180”.  (d)-(f)  Diffraction  limited  images 
recovered  by  solving  (13)  at  the  same  object  distances  as  for  (a)-(c). 


5.  DISCUSSION  AND  CONCLUSION 

We  have  utilized  a  computational  imaging  technique  along  with  a 
nonlinear  reconstruction  method  and  demonstrated  that  extended 
depth-of-field  is  possible  for  passive  millimeter  wave  imaging  even 
when  the  cubic  phase  data  is  sparsely  sampled. 

Millimeter  wave  systems  image  temperature  contrasts.  Hence, 
careful  analysis  of  noise  and  contrast  in  such  systems  in  necessary 
to  assess  the  impact  of  inserting  a  cubic  phase  element  in  the  pupil 
plane  of  the  system  and  sparsely  sampling  the  data.  More  analysis  on 


the  artifacts  introduced  by  the  cubic  phase  element  and  random  un¬ 
dersampling  in  terms  of  the  point  spread  function  will  be  discussed 

elsewhere  (see  [18]  for  more  details). 
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